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ABSTRACT 

The next generation of proposed galaxy surveys will increase the number of galaxies with photomet- 
ric redshift identifications by two orders of magnitude, drastically expanding both the redshift range 
and detection threshold from the current state of the art. Obtaining spectra for a fair sub-sample of 
this new data could be cumbersome and expensive. However, adequate calibration of the true redshift 
distribution of galaxies is vital to tapping the potential of these surveys to illuminate the processes 
of galaxy evolution, and to constrain the underlying cosmology and growth of structure. We examine 
here a promising alternative to direct spectroscopic follow up: calibration of the redshift distribution 
of photometric galaxies via cross-correlation with an overlapping spectroscopic survey whose members 
trace the same density field. We review the theory, develop a pipeline to implement the method, ap- 
ply it to mock data from N-body simulations, and examine the properties of this redshift distribution 
estimator. We demonstrate that the method is generally effective, but the estimator is weakened by 
two main factors. One is that the correlation function of the spectroscopic sample must be measured 
in many bins along the line of sight, which renders the measurement noisy and interferes with high 
quality reconstruction of the photometric redshift distribution. Also, the method is not able to dis- 
entangle the photometric redshift distribution from evolution in the bias of the photometric sample. 
We establish the impact of these factors using our mock catalogs. We conclude it may still be nec- 
essary to spectroscopically follow up a fair subsample of the photometric survey data. Nonetheless, 
it is significant that the method has been successfully implemented on mock data, and with further 
refinement it may appreciably decrease the number of spectra that will be needed to calibrate future 
surveys. 

Subject headings: galaxy surveys, weak gravitational lensing, photometric rcdshifts 



1. INTRODUCTION 

It is essential to calibrate the true redshift distribution 
of galaxies in a photometric survey if the survey is to 
be utilized to its full potential. One application of sur- 
vey data that requires a detailed understanding of the 
distribution of galaxies is weak gravitational lensing to- 
mography. The shearing of the shapes of distant galax- 
ies via weak gravitational lensing is a powerful cosmo- 
logical probe that can be used to study the distribution 
of dark matter, the nature of dark energy, the forma- 
tion of large scale structures in the universe, as well as 
fundamental properties of elementary particles and po- 
tential modifications t o the general theory of relativity 
(recen t studi es include [Hoekstra fc Ja in (2008|_ jThomas| 
et^L] ( [2009| , |Kilbinger et^ ( p0X)I| ) , [ichiki et al.| ( |2009p ). 
Cosmic shear measurements statistically examine minute 
distortions in the orientations of high redshift galaxies, 
whose shapes have been sheared by intervening dark mat- 
ter structures. Although weak gravitational lensing pro- 
vides only an integrated measure of the intervening den- 
sity field, using source populations at different redshifts 
permits some degree of three dimensional reconstruction, 
known as tomography. The distortions are small (at the 
1% level) and the intrinsic orientation of the source galax- 
ies is unknown, thus large galaxy samples are required 
to map the density field and probe the growth of density 
fluctuations with precision. Existing cosmic shear mea- 
surements have already constrained the a mplitude of the 
dark matter fluctuations at the 10% level Heymans et al. 



(2004) Rhod es et a l. (20041 Massey et al. (20051 Sem 



boloni et al. 1 2006 ) and there are many exciting galaxy 
survey proposals that will increase the available number 
of source galaxies by two orders of magnitude includ- 
ing DES, DUNE, Euclid, LSST, PanStarrs, SNAP, and 
Vista. 

Because these large galaxy surveys will have photomet- 
ric rather than spectroscopic redshift identifications, the 
community has carefully attended to fine tuning the cal- 
ibration of the photometric redshifts, minimizing biases 



and catastro p hic errors Ma ndelb aum et al. | (12008 1) , |Lim a| 
et al.| (|2008l, JFreeman et al.| pooul , |Xia etaTf ( |2009p , 
Cerdes et al.| (|2009|). Unlike experiments that use the 



galaxy positions to directly trace the underlying dark 
matter distribution, such as baryon acoustic oscillation 
studies, weak lensing analyses do not require a precise 
redshift identification for each individual source galaxy. 
It is sufficient to accurately determine the redshift dis- 
tribution of the sources. However, lensing measurements 
are extremely sens itive both to e rror and bias in the 
source distribution Huterer (20021. Attaining an accu- 
rate source distribution will be crucial if weak lensing 
measurements are to be competitive with other cosmo- 
logical probes in constraining the cosmological parame- 
ters. 

Another example where calibration of the true distri- 
bution of galaxies may be essential is in using the abun- 
dances and clustering of different galaxy populations to 
connect galaxies at late times to their potential pro- 
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genito rs at early times (as in e.g. Conroy & Wechsler 



(2009)). Such studies also utilize the luminosity tunc 
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tions of galaxies in each redshift slice, and sometimes di- 
vide these into different rest-frame color bins. To avoid 
potential systematics in inferences made about galaxy 
evolution, it will be necessary to know if some fraction 
of the population in a given photometric redshift bin is 
actually living at a different redshift, especially if there 
is an asymmetry in such errors that depends on color. 

Recently, an alternative approach to attaining an ac- 
cu rate source reds hift distribution has been proposed 
in Newman (20081. This method is similar to cross- 



correlation techni ques used in |Phillipps| ( |1985| and Mas- 
jedi et ah] ( |2006| h and the idea has also been studied 



theoretically inpchneider et al. (2006) and Bernstein fc 
Huter cr (2009). A similar technique was used in Erbcn 
et al. ~ f2009p to check the redshift distribut ion for in- 
terlop ers. Simi lar in spirit, the analysis of Quadri & 



Williams (2009) uses close angular pairs of photometric 
galaxies to constrain the photometric errors without the 
use of a spectroscopic sample. 

The cross-correlation method determines the photo- 
metric redshift distribution by utililizing the cross cor- 
relation of the galaxies in the photometric sample with 
an overlapping spectroscopic sample that traces the same 
underlying density field. One advantage of this approach 
is that the spectroscopic sample used to calibrate the 
photometric redshift distribution can be comprised of 
bright rare objects such as quasars or Luminous Red 
Galaxies (LRGs) whose spectra are relatively easy to 
obtain, and indeed may already exist in legacy data. 
Spectra could also be obtained for emission line galax- 
ies (ELGs), which are easy to follow up but may not 
represent a fair subsample. Another advantage is that 
catastrophic redshift errors in the photometry do not sys- 
tematically bias the redshift distribution estimate, they 
merely contribute to the noise. 

The cross-correlation method makes use of two ob- 
servables, the line-of-sight projected angular cross- 
correlation between the photometric and spectroscopic 
samples w ps (9), and the three dimensional autocorre- 
lation function of the spectroscopic sample £ ss (r). By 
postulating a simple proportionality between the auto- 
correlation function of the spectroscopic objects and the 
three dimensional cross-corelation function between the 
two samples £ ps (r) oc £ ss (r), it is potentially possible to 
infer a very accurate redshift distribution for the photo- 
metric sample. This assumption is guaranteed to be valid 
if the spectroscopic sample is a sub-sample of the pho- 
tometric population, but may be problematic if the two 
sets of tracers have different bias functions with respect 
to the dark matter. 

In this paper we develop a pipeline to apply the cross- 
correlation method. In section [2] we review the theory 
and explain how the method works. We highlight its 
strengths and examine potential drawbacks and system- 
atic effects. As a proof of concept, in section[3]we use the 
halo model to populate N-body simulations with mock 
photometric and spectroscopic galaxy data to quantify 
the properties of this redshift distribution estimator. We 
examine the extent to which different bias functions in- 
terfere with the reconstruction of the true distribution 
of photometric galaxies. In section [3] we discuss inherent 
tradeoffs, outline outstanding theoretical questions, and 
draw our conclusions. We leave a more detailed discus- 
sion of error propagation to the appendix. 



2. THEORY 

The goal of this paper is to demonstrate via numerical 
simulations that two spatially overlapping samples, one 
photometric and one spectroscopic, can be combined to 
infer the redshift distribution of the photometric sample 
to very high accuracy. The redshift distribution is 
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where g— ^ is the number of photometric galaxies per 
unit redshift, per steradian, and the the quantity in 
brackets is the the total number of galaxies (per stera- 
dian) in the sample, ensuring that <j) p (z) integrates to 
one. If a survey is divided into a number of redshift bins 
Zi, then 4>p( z i) gives the fraction of the total number of 
galaxies that live in the i th bin. Suppose we observe the 
angular cross-correlation function between all the photo- 
metric galaxies and the spectroscopic galaxies in a par- 
ticular bin Zi. This angular cross correlation function is 
related to the photometric redshift distribution that we 
are attempting to calibrate. 



w. 
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Here £ ps ( r ) is the three dimensional cross correlation 
function between the entire photometric sample and the 
spectroscopic galaxies that live in bin i, which is not 
observable because the redshifts of the photo-z sample 
are not known to sufficient accuracy to measure it. The 
key assumption of the cross correlation method is that 
^p S (r) tx S, ss (r), where £, ss (r), the 3D autocorrelation 
function of the spectroscopic calibrators, is observable. 
This is a reasonable assumption because on large (linear) 
scales, both £, ps (r) and £ ss (f) are related to the underly- 
ing dark matter power spectrum Af in (k) as 
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where b s and b p are the linear biases of the spectroscopic 
and photometric samples and jo is a spherical Bessel 
function. 

In asserting the proportionality, it is implicitly as- 
sumed that these biases are scale independent and that 
they evolve similarly with redshift. The former assump- 
tion is valid unless the correlation functions are being 
measured on scales smaller than ~ 1 Mpc/h, while the 
latter assumption may in fact present some difficulty un- 
less the spectroscopic objects are a fair subsample of the 
photometric population. This is because in real life, the 
photometric sample may be apparent magnitude limited. 
Thus, the population of galaxies being examined at high 
redshift may be systematically brighter, rarer, more bi- 
ased objects than those at low redshift, so the bias can 
be expected to evolve in a way that will be difficult to 
reliably calibrate. One principle goal of this paper is 
to examine the extent to which this impacts the cross 
correlation method, particularly whether the systematic 
biases involved are substantial compared to the resolu- 
tion and accuracy of the photometric distribution that is 
recovered. 
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To be very specific, on linear scales the relationship be- 
tween the cross correlation function and the (observable) 
autocorrelation function of the spectroscopic sample is 
given in equation [5] (further corrections are required for 
translinear scales). 



forms of 4> p (z) used later in this paper (see sec. 3.1 1 and 



(5) 



Thus we may write 

w ps (6,Zi)= ■^ T -r£ >ss {r{z,z i ,6))(j) p {z)dz (6) 
Jo °s(Z) 

Since b s (z) can be fit with the spectroscopic data, this 
means that we can only invert the relation to solve for 
the product b p (z)<j> p (z) in terms of observable quantities. 
This degeneracy cannot be resolved by measuring the 
angular correlation function of the photometric sample. 
We can express w pp in terms of £ ss 

w pp {6) = 

J d Zl J dz 2 M*i)M**) UP, (7) 

Changing variables to a central redshift and a Az = z\ — 
Z2, and taking note that !; ss vanishes for large Az, we 
find that 



w pp (9) oc dz 4> p {z) 
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In terms of known observables, this relation can be in- 
verted to determine the product K,(z)(j) p (z). Unfortu- 
nately, this quantity has the same direction of degeneracy 
as the product b p (z)<ft p (z). Thus the observable w pp {9) 
can only be used to improve the accuracy to which the 
product is determined, it cannot be used to break the de- 
generacy between the large scale bias and the selection 
function. In order to obtain 4> p (z) it will be necessary 
either to appeal to some model of the bias evolution or 
to find another observable that can be used to break the 
degeneracy. 

This is significant because estimators of e.g. the mean 
redshift of a sample will be affected by assuming a func- 
tional form for the bias that is incorrect . If we know we 
are recovering the product b(z)<f>(z) then in our estimator 
of the mean redshift 



^■est 



while the true z is 
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It is not unreasonable to suspect that the bias may not 
be a particularly smooth in its transitions if the sam- 
ple of galaxies accessible to a photometric survey shifts 
abruptly as some redshift threshold is crossed. One ex- 
ample of a very rapidly changing bias funct ion can be 



seen in table 2 of Padmanabhan et al. ( 2009 1 where the 
bias of a sample of LRG galaxies is computed in several 
photometric bins. At a redshift of z ~ 0.35, the bias 
jumps dramatically from 1.77 to 2.36, and back down 
again (somwhat more smoothly) to 1.9 at higher z. As 
a quick illustrative example, we take the two functional 



compute the error in z that occurs if we assume a smooth 
transition in the galaxy bias from 1.7 to 1.9 in 6 est , but 
allow the true bias &true to jump to 2.3 in between. We 
find that the fractional error in the mean redshift is 6% 
and 11% if the interval is < z < 1 and the jump is 
placed at z = 0.8. If the jump is place at z = 0.3, the 
fractional error is 0.05% and 0.2%, which is somewhat 
less significant, since there is substantially less volume 
at z = 0.3. 

3. TESTS WITH MOCK DATA 
3.1. Simulations and mock galaxy catalogs 

We use the halo model to populate N-body cold dark 
matter simulations with mock galaxies to test the cross- 
correlation method. The simulations compute the evolu- 
tion of large scale structure in a pe riodic, cubica l box of 
side 1 Gpc/h using a Tree-PM code [White) ( |2002l ). There 
are 1024 3 dark matter particles of mass o x 10 M©//i. 
The randomly generated Gaussian initial conditions are 
evolved from a starting redshift of z — 75. The Plummer 
softening is 35 kpc/h (comoving). 

Halo catalogs are constructed fro m this simulat i on us- 
ing a Friends of Friends algorithm |Davis et al.| (1985) 
with a linking length of b = 0.168 in units of the mean 
inter-particle spacing. There are approximately 7.5 mil- 
lion halos of mass greater than 5 x 10 n M©//i resolved 
with ~ 8 or more dark matter particles. The galaxy cata- 
logs are constructed by populating these halos according 
to the halo-model prescription. It is assumed that very 
small halos host no galaxies. Halos that cross a mass 
threshold M min are assumed to host one central galaxy. 
Halos of much higher mass are assumed to have formed 
through mergers of smaller halos and will host both cen- 
tral and satellite objects. The mean number of galaxies 
in a halo of mass M is given by 



<iV g ai(M)> = 0(M - M min ) 1 + 



M-M mia \ 
10M min J ( ' 

The position of the central galaxy is assumed to be at 
the halo center defined by the position of the most grav- 
itationally bound dark matter particle. The number of 
satellite galaxies in any particular halo is drawn from a 
Poisson distribution, and the satellites are assumed to 
trace the dark matter. 

We use this prescription to populate the simulation 
with several different galaxy samples. One population, 
used as the mock photometric catalog, has a relatively 
low value of M m j n = 5 x 10 11 Af© / h; these arc fairly com- 
mon galaxies with low value of the large-scale bias b p 
with respect to the dark matter. Although we know the 
true value of the rcdshifts of these objects from the sim- 
ulations, we do not use any redshift information for the 
photometric sample when testing the reconstruction al- 
gorithm developed here. The other populations we have 
created are mock spectroscopic samples, with values of 
M mln = 1 x 10 12 M e /h and 7 x 10 12 Af©//i that generate 
much rarer mock galaxies with higher biases b s with re- 
spect to the underlying dark matter. For some tests of 
the cross-correlation method, we also use a fair subsam- 
ple of 50% or 80% of the mock photometric population 
to define the spectroscopic sample. 

To test the cross-correlation method, we also need to 
impose both a selection function and an observation win- 
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Fig. 1. — The final photometric redshift distribution is affected 
by the selection function and also the survey geometry. 

dow on the photometric sample. For convenience, we 
choose to perform the reconstruction in bins of comoving 
distance Xi rather than redshift bins zj, but the method 
can be used with either choice of bins. We place the ob- 
server in the center of one face of our simulation cube, 
and assume that (s)he observes a conical volume with a 
12 degree opening angle. The cone stretches the length 
of the box (1 Gpc/h) along the line of sight. We adopt 
a toy model for the selection function, the sum of two 
Gausseans, which defines the fraction of galaxies "de- 
tected" Nk oep (x)/N(x) at each of 70 slices in comoving 
distance x through the box. 



^kecp(x) 

N(x) 
1 



"(x-X2) 2 /2o 



(12) 



We have normalized this quantity so that the maximum 
value is 1, and we refer to it as the "detection fraction" 
later in the paper. We present two models for compari- 
son in this paper, [xi, cr±, X2, ^2] = [0,0.15,0.8,0.16] and 
[0.3, 0.07, 0.7, 0.10]. In the limit of small galaxy numbers 
in the slice, it is significant that we apply the selection 
function before the geometry for reasons of cosmic vari- 
ance. The shape of the final photometric redshift distri- 
bution is affected both by the selection function and the 
geometry of of the mock observation. This is illustrated 
graphically in figure [T] which shows the first of the two 
selection functions. Once we have mock photometric and 
spectroscopic samples in place, we compute in each bin 
Xi the autocorrelation function of the spectroscopic sam- 
ple ( ss (r, and the angular cross correlation between 
the spectro-z objects in the i th bin and the entire photo- 
metric sample w v s (9,Xi)- We have used the algorithm of 



Landy and Szalay Landy & Szalay ( 1993 1 to compute the 
correlation functions. To compute w ps {9, Xi) we opt to 
measure the 3D correlation function £p S {r,Xi) an d per- 
form the integral of eqn. [(^numerically, because the mock 
observation volume is small which renders direct calcu- 
lation of w ps (9, Xi) noisy. This will not be a problem for 
surveys that cover a large fraction of the sky. 

We are eager to test how sensitive the reconstruction 
of the photometric distribution function 4> p (x) may be 
to evolution of the bias b p that is not accounted for. We 
addressed this question analytically in section [2] but it is 
useful to examine the issue with simulations to see if the 
systematic biases that occur are significant compared to 
the error bars on the re-constructed distribution. Evo- 
lution of the bias can occur because of evolution in the 
underlying large scale structure, and it can occur be- 
cause the population of photometric objects detected by 



an instrument evolves as a function of redshift, as ex- 
pected in a magnitude limited sample. Our simulation 
volume is made of a single time-slice, so we cannot exam- 
ine the effects of the first mechanism at present, but we 
expect the effects of the second mechanism to dominate 
the bias evolution. We h ave devised a simple m ethod 
inspired by the results of Vale & Ostriker (20061 and 



Conroy & Wechsler (20091 to introduce this type of bias 
evolution into our mock galaxy sample. Whereas before 
we applied the selection function in eqn [12] randomly to 
the galaxies in each slice, now we choose to keep galaxies 
preferentially that live in the largest halos. We assume 
that the brightest galaxies live in the biggest halos and 
that central galaxies are brighter than satellite galaxies 
at a given redshift. First we rank order the halos in the 
slice. We determine the number of galaxies to be kept 
from eqn. |12[ and place one in the center of each halo 
beginning with the halo of highest mass. If the number 
of galaxies exceeds the number of halos in the slice, we 
begin placing satellite galaxies. We again start with the 
largest halo and continue populating each halo with satel- 
lites until we have placed all the galaxies. This procedure 
generates a mock catalog of photometric galaxies whose 
clustering strength will depend strongly on the selection 
function. The effect of this procedure on the largescale 
bias is somewhat complicated. In the regime where only 
satellites in the lowest mass halos are being eliminated, if 
the detection fraction from eqn. 12 is high the bias will 
be close to the bias of the whole population. However 
as the detection fraction gets lower, the sample will be 
more highly biased with respect to the dark matter be- 
cause only satellites in the largest halos are being kept, 
thus weighting the largest halos more heavily than the 
smaller halos in the clustering measurement. However, 
if the detection fraction is so low that all of the satellites 
are eliminated and the population consists only of cen- 
trals, the largescale bias will be lower than for the whole 
population, because more weight from the largest halos 
has been discarded than from the smaller halos. In this 
analysis we are principally in the second regime, which 
means the large scale bias will tend to follow the shape 
of the selection function. 

This is demonstrated in fig. [2] which plots the cor- 
relation function measurement of the mock photometric 
sample in each conical slice. The top panel shows the 
result if the galaxies are eliminated randomly, and the 
bottom panel shows the results from the procedure we 
just described. The inset shows the variation in the value 
of £, PP (r) along the line of sight (where we have picked a 
particular scale, indicated with a vertical line). On the 
top we show that compared to the last slice (marked 
with blue squares), the value of £p p {r) varies by some 
tens of percent, and is a relatively flat function of the 
line-of-sight distance. The bottom inset on the other 
hand, clearly reflects the shape of the redshift distribu- 
tion function used to create the sample (plotted later in 
the bottom panel of fig. , and shows variations in the 
normalization of £, pp (r) of up to 150%. 

We emphasize that we do not expect this method to 
quantitatively capture the bias evolution in a real galaxy 
survey, but we expect it is qualitatively similar, and as 
long as the bias evolves significantly, it will be useful in 
testing the magnitude of the effect on the reconstruc- 
tion. We will refer to this method of applying the selec- 
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Fig. 2. — The 3D correlation function of the photometric sam- 
ple in each slice along the line of sight. The closest and farthest 
slice plotted are labeled with symbols. The inset shows the rela- 
tive variation among the curves at a single scale. On the top the 
there is much less variation among the curves than on the bottom. 
The inset on the bottom clearly reflects the shape of the selection 
function (and hence the bias) used to create the sample. 

tion function as the Ordered Selection Function (OSF), 
and the simple random method as the Random Selection 
Function (RSF). By comparing the reconstruction in the 
two cases, we determine how important it is to account 
for evolution in the bias. 



3.2. Pipeline 



Once the cross correlation w ps {9,Xi) and autocorrela- 
tion function ^ S a( r iXi) have been measured, a procedure 
is needed to perform the inversion of the integral of eqn. 
[6] to obt ain the line of sight distribution of the photomet- 
ric sample in some bins Xji an d to obtain error bars on 
those measurements. For a single angular scale 6 and bin 
i in spectroscopic redshift we approximate the integral in 
eqn. ?? as a discreet sum over N redshift bins (ignoring 
bias evolution for the moment). 

/•OO 

w P s{0,Xi)= / £,ss{r{x,Xi,8))<Pp{x)dx 
Jo 

N-l 

~ Xi ^V-3 ( 13 ) 

j=0 



The observation will be performed in multiple angular 
bins 9. The values of r in Xij will not be at the exact 
positions where we have made the measurement of £ ss (in 
our case in 60 bins between 5 and 50 comoving Mpc/h). 
If the r falls within the domain of our data, we interpolate 
£, ss with a spline, if it is larger than the biggest scale we 
measure, the value of £ ss is extrapolated using a fit of 
the form r~ 2 to a subset of data points (larger than 15 
Mpc/h) measured in the volume. If the value of r is less 
than the minumum value of r measured, we reject the 
entire 9 bin, since we do not trust the method for £ ss 
measured on very small scales. 

Collecting all the observed w ps values in a single vector 
w (one clement for each unique combination of Xi an d 
9), the solution for 4> p (xj) wm be obtained by inverting 
relation 

w = X-<j> (15) 

Here, X is a matrix whose elements are given by eqn. |13| 
and whose dimension is 

(#specbinsi * #thetabins) x N 

<p will be a vector of length TV, and w will be a vector 
of length (# spec bins i* # theta bins). It is significant 
that measurements at different values of 9 can mix in the 
inversion; since they are correlated it is important that 
they do so. We do not know the pre-factor b p /b s , but 
as long as it is assumed to be constant (as in the special 
case where the spectroscopic sample is a fair subsample 
of the photometric population) this is not relevant. The 
reconstructed 4> P (xj) wm not be correctly normalized af- 
ter the matrix inversion because there is no information 



in eqn. 13 about the total number of galaxies in the pho- 



tometric sample. The normalization of 4> P {Xj) wm be set 
by requiring that it integrate to L. 

In principle, solving equation [15] for (f> p (xj) is a sim- 
ple matter of inverting a matrix, but in practice doing 
so is numerically unstable and furthermore errors in the 
observables w and X must be accounted for in the solu- 
tion. Since to is a projection through the same dark mat- 
ter structures traced by X, there will be non-negligible 
correlations between the errors in the two observables. 
To extract 4> p (x)> we write down the expression for the 
X 2 statistic, and minimize it. Typically we are attempt- 
ing to reconstruct <p p (xj) m as many bins N as possible. 
There is often insufficient information in the correlation 
functions to completely constrain <j) p in the chosen num- 
ber of bins, thus it becomes necessary to stabilize against 
solutions that have highly anti-corrclatcd adjacent bins. 
Since it is reasonable to assume that the true distribu- 
tion will not be highly oscillatory, we adopt a smoothness 
prior to regularize the solution, and add it to our expres- 
sion for x 2 below. 

X 2 = (w - X<f>) T C^(w - X0) + A 4> T B T B<p (16) 

Here C<£ is the covariance matrix incorporating errors in 
both w and X, and the subscript denotes that it is an ex- 
plicit function of <p. For purposes of the present analysis, 
however, we will assume that the spectroscopic correla- 
tion function is perfectly determined, and we will prop- 
agate errors from w ps only, specifically CT 1 = C _1 = 



Xij — b 



(Oxi) 2 + (Xj~X*f) (14) 



inverse covariance matrix of 



(9,Xi)- This is not 



good approximation as we will soon show, and we refer 
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the reader to appendix [A] for a description of the full 
error propagation from C^ 1 into errors on <p p {xj)- The 
reader should consider all errors reported on <ft p (xj) as 
lower limits on the total error. 

The regularizatio n scheme w e have adopted is de- 
scribed in detail in Press (20021 section 18.5. The inten- 
tion is to add a term to x 2 that gets large when neigh- 
boring points have widely different values. Minimization 
of x 2 w ih then tend to solutions that do not have anti- 
correlated neighboring points. The matrix B is given 
by 



B 



/-l 1 00 -0\ 
-110 •••0 



•••00 0-1 1 
V •••00 -1 1/ 



(17) 



Note that B has one fewer row than column. The fac- 
tor A should be chosen such that the first and second 
terms in x 2 contribute roughly equal weight. This can 
be approximately arranged if A is taken to be 



A = 



Tr(X T C~ 1 X) 
Tr(B T B) 



(18) 



but in practice, the weight in x 2 will also depend on the 
solution and the values of w ps . We find that the quality of 
the reconstruction depends on how well the two terms in 
eqn. 16 are balanced. Since this depends on the answer, 



we opt to refine the value of lambda and re-compute the 
solution iteratively as follows. 

• Compute a tolerance parameter defined from the 
two terms in x 2 &s tol=l-(term 1/term 2). 

• If tol > 0, A is increased by a factor of 10. If tol 
< we decrease it by 10. 

• We recompute the solution and the tol 

• If the absolute value of tol has decreased, we repeat 
the refinement of lambda, otherwise we exit and 
keep the previous value of the solution 

In practice, the procedure usually requires only 1 refine- 
ment of A. We observe that changing the algorithm to 
have smaller steps in lambda (e.g. a factor of 2 rather 
than 10) does not improve the solution, and occasionally 
over-smooths it. We emphasize that while we have iden- 
tified an algorithm that works, we have not optimized 
the application of a smoothness prior. Since the solution 
is moderately sensitive to the smoothing, care should be 
taken to understand the properties of the smoothing be- 
fore applying this technique to real data. We have not 
studied the effects of different smoothing algorithms be- 
cause it is likely that the need for a smoothness prior will 
be eliminated in future analyses by using the photomet- 
ric probability distributions as a prior instead. Since we 
have no mock photometric redshift probabilities, exam- 
ining that technique is beyond the scope of this paper, 
but will be the subject of further research. 

To minimize x 2 we take the derivative with respect to 
4> and set it equal to zero. Taking C _1 to be independent 
of 4> and noting that it is symmetric we find 

- 2w T C- 1 X + 2<j) T X T C- 1 X + 2\<j) T B T B = (19) 



X T C- X w = [X T C- X X + XB T B] </> (20) 
= [X T C- l X + \B T B] _1 X T C- 1 w (21) 
The covariance matrix of the recovered cb is given by 



Cov[0] oc [X T C _1 X + XB T B] 



(22) 



To recover the constant of proportionality, we will need 
to rescale the elements of the covariance matrix by the 
square of the factor used to rescale 4> (since we renor- 
malize such that cj)(x) integrates to 1). We caution that 
all matrix multiplications above are finite sum approxi- 
mations to integral quantities, so care must be taken to 
ensure that phi and its error bars come out to scale. For 
clarity, let us refer to w ps (zi) as a function of a continu- 
ous variable Zi (which will label the spectroscopic bins). 
£,(zi,Xi) will be a function of both Zi and another con- 
tinuous variable x% (which will label the bins in which <f> 
is reconstructed) . We are considering a single value of 6, 
and z and x can represent either redshifts or comoving 
distances, we simply use both letters to distinguish them. 
Ignoring regularization, the continuous expression for x 2 
is then 



X 



dz\dz2 



U'v 



Ol) - / dxKt>(Xi)£(zuXi) 



C 



' 1 (Z1,Z 2 ) (w ps 



O2) - / cfe<HX2)£(Z2,X2) 



(23) 



To minimize x 2 we must set the functional derivative 
Sx 2 /S4>(x) — 0- We leave the details to the reader, but 
note that in eqn [19} the factors of Az that correspond to 
integrals over z cancel but the factor of Ax does not. 

3.3. Redshift Distribution Reconstruction 

In this section we present a series of reconstructions 
that examine various aspects of the problem and iden- 
tify the potential difficulties in applying this method. 
We begin in fig. [3] with the simplest case and gradually 
add complexity. The heavy solid line shows the theo- 
retical redshift distribution (selection function + geome- 
try) that has been applied mock photometric catalog of 

This redshift 
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galaxies with M m i n = 5xlO n in eqn. 
distribution corresponds to the selection function shown 
on the left of fig. [I] We have chosen g alax ies randomly 
(the RSF method described in section 3.1l in applying 
the selection function. The resulting catalog has 54,000 
galaxies. We have measured the detection fraction in 
each slice for this particular realization, and plot it with 
a thin wavy line in fig. [3] The spectroscopic calibrating 
sample is taken to be a different population of galaxies 
with M m i n = 7x 10 12 yeilding around 5600 galaxies, i.e. 
a highly biased tracer of the density field compared to 
the photometric sample. We assume both the observ- 
ables are measured with perfect accuracy. Since we have 
assumed perfect knowledge of the spectroscopic correla- 
tion function, we opt to measure it in the entire light 
cone, and us the result in each of the i spectroscopic bins 
£, ss (r,Xi) = £whole volume {r)- This is only justified because 
there is no evolution in the light cone, and the calibrators 
are not affected by the selection function which changes 
along the line of sight. We perform the reconstruction by 
computing the matrices -X", w, C , and B and using 
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aai Theory 
— Realization 
■ ■ Reconstruction 




0.2 0.4 0.6 0.8 

Comoving line of sight distance \ (Gpc/h) 



Fig. 3. — The photometric distribution </> p (x) an d its reconstruc- 
tion. The x axis is comoving line of sight form the observer, who 
is located at the origin. The thick solid line shows the theoretical 
value of the distribution, the thin wavy line shows the particu- 
lar realization in this simulation volume, and the points show the 
reconstructed solution. 



them in eqn. [2l] The result is plotted in fig. [3] with 
points. Notice that the reconstruction follows the values 
of this realization rather than the theoretical redshift dis- 
tribution used to create the sample. 

Fig. [4] reveals a more realistic picture. The lines and 
points in fig. [4] are all the same as in fig. [3] We now 
measure the autocorrelation of the spectroscopic sam- 
ple £ S s(r, Xi) m each of the conic sections, and use those 
measurements to perform the reconstruction. The top 
and bottom panels show two different choices of selection 
function, whose parameters are marked in the caption. 
The top is difficult to reconstruct because it peaks at 
X = where the is no volume in the mock observation, 
and the bottom is a challenge because the feature coin- 
cides with the bin spacing, which will be a problem to 
reconstruct because of the smoothness prior. We have 
drastically decreased the number of bins, so that there 
are a reasonable number of calibrators in most of the 
bins. The normalization criterion, which comes from the 
condition that <f) p integrate to 1, becomes more sensitive 
to noise in the reconstruction when the number of bins 
is decreased. In this plot and the plots that follow we set 
the normalization by hand so that we can illustrate other 
points; we let the maximum value of the reconstruction 
(points) equal to maximum value of the theoretical in- 
put 4> (thick solid line). While the reconstructions in 
fig. [4] show the correct trends, they are not sufficiently 
high quality to recover the bimodal behavior in the re- 
constructed selection function. 

We now bri efly discuss the error bars in fig. [4]obtained 
from eqn. 22 The Cov[</>] matrix is not diagonal; we con- 
servatively report 1/ -y/cov -1 ^]^ as the error on the i th 
bin. These error bars represent a lower bound on the er- 



ror because we have only propagated error from w pa and 



ing the observation into a number of bins and bootstrap- 
ping. Here, however, we opt to follow Peebles (1980) and 
assume that on these scales the errors are dominated by 
Poisson noise. In this case, C _1 is diagonal, and its ele- 
ments are given by 

1 



5 £ p X«Ax s nsin(0)A0 (24) 



pairs 



We use survey parameters appropriate to large upcoming 
missions. We take h s = 1 X 10 5 galaxies per (Gpc/h) 3 , 
S p = 100 photometric galaxies per square arcmin, and 
ft = 20, 000 square degrees on the sky. Ax s is the width 
of the bin of spectroscopic data. The resulting error bars 
in [4] and those that follow display a surprising trend. 
Even though the number of pairs is dramatically larger 
for bins at farther comoving distance, the reconstruc- 
tion of the photometric distribution is less accurate in 
the most distant bins. The key to understanding this 
trend is that for a fixed angular scale 9, the physical 
scale being probed in distant bins is much larger than 
at nearby distances. Since the correlation function is a 
steeply dropping function of separation, the impact of 
shrinking values in X of eqn. [22] dominates over the 
growing number of pairs in C -1 . By measuring w ps on 
smaller angular scales the errors in the farthest bins can 
be improved, but there is a limit to how far this can 
be pushed because measurements at small angular scales 
will send the near field into the trans-linear and 1-halo 
regime, for which there are corrections to the underlying 
assumption that £ ps oc £ ss . 

In figure [5] we switch to a larger calibration sample with 
M m i n = 1 x 10 12 . We see that the reconstruction now 
captures the bimodal behavior, but both the resolution 
(bin spacing) and errors are modest and the smoothing 
is still evident in the last bin. Notice that the error bars 
have increased: this is because this calibration sample 
has lower bias, so the elem ents in X in the covariance 
matrix of <f>(x) (eqn. 22 1 are all smaller. This is an 



not from £ ss . The matrix C in eqn. 22 is the inverse 
covariance matrix of the cross-correlation measurement 
iu, and with sufficient volume can be estimated by divid- 



illusion however. Had we propagated the error from £ ss , 
it would contribute larger errors to figure H than to [5] 
because the mcasurmcnt is noisier for the smaller sample. 
In moving to the larger sample, the number of spectra 
required has increased by an order of magnitude, and is 
now the same size as the mock photometric catalog (after 
the selection function has been applied). In summary, 
comparison of fig. [3j fig. [4] and fig. [5] demonstrates that 
£ss(x) must be reasonably well determined for the shape 
of the reconstructed photometric distribution to be well 
captured. It is worth noting that we have used no priors 
in the determination of £ S s(x), improving and applying 
theoretical priors may yield significantly better results 
with fewer spectra. 

The fact that the more numerous calibrators generated 
such improved results suggested to us that perhaps the 
reconstruction with the rarer calibrators could be im- 
proved by simply fitting the spectroscopic observations 
with a 2-parameter power law, and performing the re- 
construction in that way. This appears to discard too 
much information contained in £ ss (x)- The result was 
visually similar to fig. [4j namely, the bimodal behavior 
in the distribution is lost. 

We add another layer of complexity in fig. [6] We ap- 
ply our selection function in such a way that the bias of 
the remaining objects will evolve along the line of sight 



8 



A.E. Schulz 



0.02 



^ Theory 
— Realization 
□ o Reconstruclion 




0.2 0.4 0.S 0.8 

Comoving line of sight distance \ (Gpc/h) 



0.02 



^ Theory 
— Realization 
m a Reconstruction 




0.2 0.4 0.S 0.8 

Comoving line of sight distance \ (Gpc/h) 




0.2 0.4 0.S 0.8 

Comoving line of sight distance \ (Gpc/h) 
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Fig. 4. — The photometric distribution <j> v (x) an d its re- 
construction. The top and bottom panels show two differ- 
ent selection function choices. Each is the sum of two Gaus- 
sians with [xi, 01 , X2, ""2] = [0.0,0.15,0.8,0.16] on the top and 
[0.3,0.07,0.7,0.10] on the bottom. The errors come from Pois- 
son error in the cross-correlation measurement. The spectrosc opic 
sample is made of rare objects with M m ; n = 7 X 10 12 in cqn. |ll| 
There are very few objects in each conic section so the correla- 
tion functions are poorly measured, and this reconstruction is not 
capturing the bimodal behavior. 



(the OSF method described in section 3.1 1. We continue 



using the less rare calibrating sample with halo model 
A/min = 1 x 10 12 for this illustrative proof of concept. Al- 
though it is not very statistically significant, the eye can 
pick out a trend in the reconstruction: the reconstructed 
distribution is suppressed in areas of low detection frac- 
tion. Since the bias is tracking the detection fraction, 
the areas of low detection are enhanced less than the ar- 
eas of high detection, thus they appear suppressed when 
we normalize to the highest point. At this level of accu- 
racy, it is unlikely that this systematic will dominate the 
error in tomographic analyses, however for much larger 
surveys it could be a significant concern. 

In fig. 7] we show how the situation is altered if a fair 
subsample of the photometric population is used instead 
of a rarer biased tracer population. This is a special 
case because the bias of the calibrators will evolve with 
redshift identically to the photometric sample. We show 
two different subsamples in this plot, 50% and 80% of the 
photometric catalog, and we test the reconstruction for 
the distribution on the right in figs. 4|5|6 We see that for 



Fig. 5. — The photometric distribution (pp(x) an d its reconstruc- 
tion. Now the spectroscopic population has Af m i n = 1 X 10 12 and 
is roughly the same size as the photometric one though it is more 
biased and comprised of different galaxies. Because the correlation 
function has been measured well, the reconstruction recaptures the 
bimodal behavior. 

a survey of this volume, 50% is too small to capture the 
bimodal behavior in the reconstruction, but with 80%, 
the reconstruction works well. Indeed the systematic off- 
set in figure [6] is remedied and the reconstruction follows 
the realization to much better accuracy. For larger sur- 
veys the fraction needed will be smaller than indicated 
here, because they will be able to measure the spectro- 
scopic correlation functions with greater accuracy. For 
surveys large enough not to be limited by noisy spec- 
troscopic correlation functions, it may be necessary to 
calibrate with a fair subsample of galaxies to avoid sys- 
tematic bias in the redshift distribution that comes from 
evolution in the bias of the photometric sample. 

4. DISCUSSION 

We have outlined the theory behind the cross- 
correlation method for calibrating the redshift distribu- 
tion of objects with photometric redshifts, and developed 
a pipeline than can be used to apply the method to sur- 
vey data. We have created mock simulations to test the 
pipeline. We have succeeded in reconstructing the red- 
shift distribution of the mock photometric galaxies us- 
ing the angular cross correlation of these galaxies with 
an overlapping spectroscopic sample (whose redshifts are 
known). We have not used any redshift information 
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Fig. 6. — The photometric distribution p,,(x) an d its recon- 
struction. Here the selection function applied to the photometric 
sample causes the bias to evolve as a function of comoving distance 
from the obsever. The calibrators are the M m i n = 1 X 10 12 sam- 
ple, whose bias does not evolve. The result is a systematic shift of 
marginal significance in areas of low detection fraction. 



about the photometric sample. We have demonstrated 
the validity of the method. We have also identified the 
aspects that are likely to be the limiting factors in relying 
upon this method to provide accurate redshift distribu- 
tion information. These limiting factors are 1) that the 
spectroscopic sample must be binned along the line of 
sight, causing their correlation functions to be noisy and 
interfering with the reconstruction, and 2) that the bias 
evolution of the photometric sample cannot be disentan- 
gled from the redshift distribution that is reconstructed, 
which may force the necessity for the follow up of a fair 
subsample of galaxies. Improved modeling of theoretical 
priors could yield large dividends if these factors can be 
mitigated. 

The analysis has revealed a number of trade-offs that 
exist in application of this method. For a given set of 
calibrators, there is a trade-off between the resolution 
(bin spacing) of the redshift distribution reconstruction, 
and the error bars on any individual point. Thus if a 
population of catastrophic outlier were discovered, for 
example, it would be interesting to investigate whether 
it is more important to know how many there are, or at 
exactly which redshift they lie. We also find a trade-off 
between the number of calibrators and the quality of the 



Fig. 7. — A fair subsample is used to reconstruct <p p (x) instead 



of the M„ 



1 X 10 sample. The bias of the subsample evolves 



identically to the bias o the photometric sample. The reconstruc- 
tion is studied only for the distribution on the right of figs. [4] [5] 
and [6] The top panel shows a subsample of 50%, which is not suf- 
ficient to capture the bimodal behavior, and 80% on the bottom, 
which does well, and shows no systematic offset due to evolving 
bias. 



reconstruction. As the number of available spectra are 
decreased, the spectroscopic redshift bins need to widen 
to maintain equivalent signal strength. This in turn af- 
fects how finely the reconstructed redshift distribution 
can be sampled. Bimodal behavior may be lost, and 
with insufficient priors (such as the smoothness prior we 
implemented here), false bimodal behavior may appear 
in the form of anti-correlated adjacent points. 

Widening the bins to compensate for fewer spectra 
also introduces another difficulty. Recall that the er- 
ror in the angular cross-correlation function goes as 



Sw 



ps 



'pairs ■ 



This means that 



Sw 



ps 



( ps 



pairs 



-i -1 



(25) 



For the purposes of this order of magnitude argument, 
suppose 4>p were constant, then to integrate to 1 requires 
<^p(x) l/AXrb- The subscript rb is used to indicate the 
width of the reconstruction bin (not the spectroscopic 
bins). When we remove it from the integral we are left 
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with an expression that is w p (r p ) 



Sw„ 



ID 



ps 



pans 



~] -1 



AXrb 



(26) 



pairs 



If the reconstruction bin is widened by a factor of two, the 
number density of photometric galaxies will have to be 
increased by a factor of 4 to preserve the same accuracy 
in the cross correlation measurement. Therefore in this 
method there is also a trade-off between the number of 
spectra that the survey can afford versus the number of 
photometric galaxies they can afford. 

In this analysis we have not made any use of the pho- 
tometric redshifts, which although not sufficiently accu- 
rate to determine the true redshift distribution of the 
population, still contain signifi cant information abo ut it. 
Modern techniques such as in |Gerdes et a l. (2009) have 
made it possible to assign a probability distribution for 
the redshift of each individual galaxy in a photometric 
survey, rather than a single best estimate and error bar. 
We propose that combining these probabilities for all the 
galaxies in a given redshift bin constitutes a reasonably 
powerful prior, that can take the place of the smoothing 
we have introduced in this analysis. This is fortunate be- 
cause the solution is frequently wrecked by the smooth- 
ing criterion, although the analysis cannot be performed 
without it. Another detail that we leave for a future 
study is that the spectroscopic sample need not be com- 
prised of a single rare population, it is quite conceivable 
to target certain regions of redshift space that are known 
to be problematic more heavily. It must also be possible 
to use less rare tracers in regions with smaller volume. 
We leave such optimization to the future. 

There are a few important factors that we have ne- 
glected in this a nalysis. As pointed out by Bernstein 



& Huterer (20091, weak gravitational lcnsing will induce 
correlations between the positions of calibrator galaxies 
in the foreground with photometric galaxies in the back- 
ground, and vice versa. This will need to be carefully 



controlled for the method to reliably calibrate redshift 
distributions. We also have not mentioned the complica- 
tion of the integr al constraint, which as discussed in e.g. 
Huff et al. ( 2007 ) can lead to very significant errors when 
the volume and the scales in the correlation function are 
of comparable size. This may be as significant an issue as 
evolution in the large scale bias of the photometric pop- 
ulation, though it may be mitigated if integral constraint 
errors are correlated between photometric and spectro- 
sco pic samples. Fortunately, i mproved estimators exist 
(in Padmanabhan et al.| ( |2007[ ) for example) and should 
certainly be incorporated into the pipeline. 

Having now demonstrated that the method is viable, 
with further refinement the cross-correlation method ap- 
plied in conjunction with direct follow up surveys may 
significantly reduce the number of spectra that are re- 
quired to calibrate photometric redshift distributions to 
the desired accuracy. There are a few other advantages 
that we have not touched upon in detail. As we have 
shown, the cross correlation method can be used to cali- 
brate the redshift distribution all the way along the line 
of sight, and as such is uniquely suited to detection and 
calibration of catastrophic redshift errors, even if these 
errors are so rare that they are missed by conventional 
follow up. Also, the reconstruction should not be ad- 
versely affected by redshift deserts, regions where no 
spectroscopic redshifts are available. We are optimistic 
that this technique will provide a useful complementary 
approach to conventional calibration techniques, and ev- 
ery effort should be made to refine the method further. 
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APPENDIX 
ERROR PROPAGATION 

The treatment in this paper proceeds by assuming that the correlation function of the spectroscopic calibrating 
sample is perfectly determined. However we show that since the tracers are rare, the measurements of £ ss can be quite 
noisy due to binning finely along the line if sight. Error in the measurement of £ ss may ultimately dominate the error 
budget, and it is therefore important to lay out the procedure for properly incorporating it. Ignoring the regularization 
term, the expression for x 2 is 

X 2 = (w-X<fr) T C^(w-Xcj>) (Al) 
Ccf, is the covariance matrix that includes errors in both w and X. The elements of the covariance matrix are 

C^ij = C y + ]T C' ijk cf> k + <P k C"vki^ 1 (A2) 

k k,l 

which depend explicitly on the solution <fi. This renders the problem non-linear and will require iterative numerical 
methods to minimize x 2 ■ To cut down on the proliferation of indices, the expressions in C are written out below for 
the case of a single angular bin 0. 

Cij = (e w {xi)e w {Xj)) (A3) 

C' ijk = (e w {Xi)et{Xj,Xk)) (A4) 

C'ijki = (edXi,Xj)et{Xk,Xl)) (A5) 

To generalize to multiple theta bins, let the index i in e.g. e w (xi) run over each unique combination of 9x, rather than 
just over %. In the expressions for e^(xi, Xj) it is worth mentioning that in general i and j run over a different number 
of bins, and k and / run over the number of reconstruction bins N. 
It is often useful in numerical minimization to have an expression for the derivative of % 2 . This is 

^ =2(w - X^C-H-X) + {w - X^C- 1 f C" 1 ^ - X0) (A6) 

To simplify the expression we have suppressed the indicies, but note that the derivative of C is a rank 3 object, with 
two dimensions of the length of w, and one dimension of the length of <j> (i.e. N). Since C is a symmetric matrix, the 
derivative is given by 

^^ + 2^^, (A7) 



